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The non-Markovian stochastic dynamics involving Levy flights and a potential in the form of 
a harmonic and non-linear oscillator is discussed. The subordination technique is applied and 
the memory effects, which are nonhomogeneous, are taken into account by a position-dependent 
subordinator. In the non-linear case, the asymptotic stationary states are found. The relaxation 
pattern to the stationary state is derived for the quadratic potential: the density decays like a linear 
combination of the Mittag-Leffler functions. It is demonstrated that in the latter case the density 
distribution satisfies a fractional Fokker-Planck equation. The densities for the non-linear oscillator 
reveal a complex picture, qualitatively dependent on the potential strength, and the relaxation 
pattern is exponential at large time. 


I. INTRODUCTION 


The presence of long distribution tails in stochastic systems means a qualitatively different picture than a familiar 
Gaussian process. In particular, the variance of the process value x diverges and the ordinary central limit theorem 
does not apply. Instead of the normal distribution, the power-law form of the asymptotic density distribution, 

(0 < a < 2), emerges. Processes characterised by such tails are stable and called Levy flights. They are typical for 
complex systems and frequently observed in many areas of science Bi- Processes involving Levy flights are well 
suited to describe phenomena far from equilibrium. Dynamics of a particle subjected to the Levy stable noise can be 
described by a Langevin equation. 


d 

dx{t) = ——V{x)dt -\- T]{dt), 


( 1 ) 


where V{x) is a potential and the increments r]{dt) obey the Levy stable statistics with the stability index a. Solution 
of the corresponding Fokker-Planck equation is simple for the linear case - when V{x) =const, oc a; or oc - and the 
resulting density reveals the same statistics as the driving noise 0. That equation, in the presence of the Levy flights 
determined by a stability index a, contains a fractional derivative defined by a Fourier transform [^, 
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For the quadratic potential, the system converges with time to a stationary state. The stationary solutions in the 
asymptotic limit have been also found for strongly nonlinear ^sterns, V{x) (x x'^ (7 > 2), and the steep slope of the 
potential modifies the distribution making the variance finite 5,@. This analysis was generalised to a two-dimensional 
problem Q and to systems driven by a multiplicative noise 8 |; the spectral analysis was performed for transport in 
an inhomogeneous medium Q. 

The above approach is Markovian (Eq.ffl is local in time) and does not take into account memory effects emerging 
in complex [lOl| and disordered systems |ll| where, as a result of traps and faults, the particle experiences periods 
of long rests. Then the process can be formalised by a jump-size and waiting-time distribution in the framework of 
the continuous time random walk (CTRW) [l^ . The problem without any potential is well known: the variance - 
existing for the Gaussian case - rises slower than linearly, which means a subdiffusion, and the characteristic function 
of the density distribution decays according to the Mittag-Leffler pattern [l^. The diffusive behaviour for which 
the variance rises with time slower or faster than linearly is often called an anomalous transport, in contrast to a 
normal diffusion. Processes characterised by the anomalous transport - e.g. the flow in porous media where fractional 
space derivatives model large motions through highly conductive layers or fractures while fractional time derivatives 
describe motionless particles [l3 | - can be conveniently handled by the subordination procedure [iM3- In that 
description, the Markovian dynamics proceeds in an operational time and is given by Eq. o (a subordinated process), 
whereas the memory effects are modelled by another Langevin equation for the physical time which is a non-decreasing 
random variable, determined by a one-sided density distribution ^ (a subordinator). In the presence of a potential, 
the process o has a stationary limit and the subordinator determines a relaxation to that stationary state. Long 
tails of the waiting-time distribution mean a deviation from the exponential Debye relaxation. Non-exponential forms 
of the relaxation are common in complex systems: viscoelastic materials, synthetic polymers and biopolymers [l8l| . 
They emerge in supercooled liquids near the glass transition temperature and in amorphous polymers [l9|; a typical 
decay form of the relaxation processes is a Mittag-Leffler function (Cole-Cole) [l^. Diversity of systems exhibiting 
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the non-exponential relaxation traces back to the universality present in complex systems: the empirical relaxation 
laws reflect a behaviour which is independent of the details of examined systems [^. From that point of view, the 
non-exponential relaxation can be understood as a result of a strong interaction between diverse elementary units. 

In a nonhomogeneous medium, one can expect that the waiting time depends on the position. The corresponding 
CTRW description is coupled [l^ ; its Markovian version contains the waiting-time distribution in a Poissonian form, 
characterised by a variable jumping rate , and resolves itself to a Fokker-Planck equation with a multiplicative noise 
[^ . In the disordered systems, many forms of the disorder may emerge [TTI | and the dynamics becomes complicated if 
one assumes that the trap positions slowly evolve with time (a quenched disorder), then the successive trapping times 
are correlated. The position-dependence of memory effects can be taken into account by applying a subordination 
technique and such a model has been recently proposed [l^. One can expect that a large density of traps or regular 
structures in the phase space the trajectories stick to invoke a long resting time at specific regions. The subordination 
approach takes into account that position-dependence by introducing a variable intensity of the random generator 
of the physical time, given by a non-negative function g{x). Then Eq.(IT]) describes the subordinated process and the 
dynamics is defined by the following set of equations, 

dx(r) = Fd(x)dT + rj^dr) 

dt{T) = g{x)i{dT), (3) 

where the increments g{dt) and ^{dt) obey the Levy stable statistics with the stability indexes 0 < a < 2 and 
0 < /3 < 1, respectively. The process ^(t) is non-decreasing and biased by the medium properties: if, in particular, 
g{x) is large, the particle stays at x for a relatively long time, compared to the region where g{x) is small. Eg. ([5]) 
can be approximated by an ordinary subordination process if we decouple the random and the position-dependent 
ingredient; this can be accomplished by a two-step procedure [l^. First, let us assume for the moment that ^ is 
given by a one-sided density with finite moments (e.g. an exponential), approximate ^(dr) by the average, {^)dT, 
and evaluate the operational time increment, Ar = {{^)g{x))~^At. Since the first equation ([3|) can be discretized 
as Ax = Fci{x)At + the system ([3]), after inserting Ar, resolves itself to the Langevin equation with a 

multiplicative noise in the ltd interpretation which defines a Markovian process. The above equation incorporates the 
spatial distribution of traps but neglects the memory effects which are important if /3 < 1. Then, as a second step, we 
take into account those effects subordinating the Markovian process to the random time determined by /3 < 1. The 
final set of equations reads, 

dx{T) = Fd{x)v{x)dT + i/(x)^/“? 7 (dr) 

dt{T) = ^{dr), ( 4 ) 

where i/(x) = [(05(2^)]”^! assume in the following {^) = 1. The process described by Eq.(|3]) was studied for Fd = 0 
[^ : it corresponds to a Markovian CTRW with a variable jumping rate, subordinated to a random time, and resolves 
itself to a fractional Fokker-Planck equation with a variable diffusion coefficient. The system (|3]) is used, for example, 
to describe such transport phenomena as advection, molecular diffusion and hydrodynamic dispersion, where the first 
equation (jd]) is related to the advection-dispersion equation whereas the subordination accounts for memory [2^ . 

In the presence of an external potential, Eq. © predicts a convergence to a stationary density distribution. However, 
the stationary state is difficult to derive and known in its exact form only for few potentials (3 • All the more difficult 
is a problem of the relaxation to that state and the time-dependent solutions are unknown even for the Markovian 
case. It is unclear, in particular, whether the Mittag-Leffler pattern is still valid when a process with a nonlinear 
deterministic force is subordinated to the random time. In the present paper, we address that relaxation problem and 
consider a more general case of the position-dependent memory. We find expressions for the asymptotic stationary 
density for the potential in a power-law form and derive the density for the subordinated process, as a function of 
the operational time, for a quadratic potential (Sec.II). Sec.Ill is devoted to the time-dependent solutions of Eq.dS]), 
Q where the relaxation function is derived and the influence of the nonhomogeneity effects on the decay rate of the 
density is discussed. 


II. STATIONARY STATES AND THE SUBORDINATED PROCESS 

We consider a stochastic system under the influence of an non-linear deterministic force, 

Fd = -A|xPsign(x), (5) 

where A > 0 and 7 =const. This form of the stimulation is a natural generalisation of an ordinary harmonic oscillator 
and often discussed in connection with the chaotic dynamics [ 2 ^ . Hence the first equation (|3]) takes the form, 

dx(r) = —X\xp i'{x)sign{x)dT + iy{x)^^°‘f]{dT), 


( 6 ) 
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and corresponds to a Fokker-Planck equation 


la 


^ = A|:[|x|V(.)sig„(.)p„(x..)l + 

where the fractional derivative is defined by Eq.([2]). The equation for the characteristic function reads 


( 7 ) 


—Po(fc, r) = -Ak—J'Hxp~^iy(x)po(x, r)] - \k\°‘T[v{x)po{x, r)]. ( 8 ) 

Eq.® cannot be analytically solved, in general, and presence of the multiplicative factor ly^x) in the last term poses 
the main difficulty. We restrict our considerations to the asymptotic regime in positions, which corresponds to small 
wave numbers, and derive a stationary solution. Due to the factor |fc|“ in the last term of Eq.®, the leading term in 
the solution is of the order |fc|“ implying the asymptotic solution in a scaling form, 

Poix,T) = a{T)po{a{T)x) = a(T)"“ (9) 


where a =const; the above solution is valid if 0 < a < 2. Moreover, it is sufficient to keep only the lowest term, k^, 
in the fc—expansion of the Fourier transform argument and the transform resolves itself to the integral 


iy{x)po{x, T)dx. 


( 10 ) 


For the evaluation of the argument of the first Fourier transform an additional assumption is necessary: we require 
that the asymptotic form of g{x) is algebraic, g{x) = \x\^ , and that argument reads a{T)~°‘\x\^~^~°‘~'^ for large \x\. 
The last term in Eq.([ 8 ]), in turn, follows from the integral (ITOll . 

- \k\°‘^[v{x)po{x,T)] = -\k\‘^{v{x)){T), (11) 

and this expression is valid for any time. Now, we take the limit of large time and look for the stationary solution: 
Pq(x, t —>■ oo) = (a:) and a(r —>■ oo) = Og. In that limit, the first transform on the rhs of Eq.® reads 

f{k)=N-^F[Na:^\x\-^-\ (12) 

where N = {\x\^~^v{x))~^ and 5 = —7 + 0 + d + 1 (0 < <5 < 2). The transform argument is normalised to u nity 
and the transform itself can be simply evaluated up to the first order in k by means of the following lemma [28|. 
If W{x) is a normalised and symmetric density distribution and dW{x') ^ ba~^x~°‘ for a; —>■ 00 (6 > 0), then 
1 — W{k) ~ 6 c|fc|“ for A: —0, where 


TT 

r(a + 1) sin(a 7 r/ 2 ) 


(13) 


Application of the above lemma to the present problem yields f{k) = N~^ — ca~°‘\k\^ and inserting of the evaluated 
expressions to Eq.® produces Fourier transform from the stationary density which contains terms corresponding to 
\k\°‘ and The equation is satisfied for 5 = a and the parameter a is determined: d = a + 7 — 1, where 7 = 7 — 0. 
Moreover, we obtain 




Att 


1/a 


(iA(a:))r(a) sin(Q!7r/2) 


and the final asymptotic solution follows from Eq.®, 


(14) 


p\^\x) = ^ r(a) sin(a 7 r/ 2 )|a:| “ (15) 

Att 

The above solution exists if 7 satisfies some conditions: it is normalisable for 7 > 1 — a + 0 and {v{x)) is finite for 
7 > 1 — a. The x—dependence of pg‘^\x) is exact for any d but {v{x)) depends on the density at small |x| and cannot 
be determined without additional assumptions. For 0 = 0, Eq. m is exact in a wide range of 7 , including negative 
values, and the condition 7 > 1 — a is required to get a stationary solution [^ . 
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Our aim is to establish how the system under the force field ([5]) approaches the stationary state (flSl) . i.e. we have 
to derive the time-dependence of the density Po{x,t). Its asymptotic form can be analytically derived for a linear 
case. From now on, we restrict our analysis to a power-law form of g{x), 

g{x) = |a;|® (6» > -a), (16) 

which is well suited to describe diffusion on fractals [30l| . It is encountered in geology where a power-law distribution 
of fracture lengths is responsible for transport in a rock ; this self-similar structure of a fracture and fault network 
is characterised by a fractal dimension and determines a pattern of water and steam penetration in the rock . The 
parameter 9 estimates a degree of the memory nonhomogeneity: for 0 < 0, g{x) is small near the origin which means 
long waiting times whereas for 0 > 0 transport is inhibited at large distances. The deterministic force is linear when 
7 = 14-0 and Eq. 0 takes the form, 


dpo{x,T) 

dr 


= + 


d°^[\x\ %(a:,T)] 


(17) 


that leads, after the Fourier transformation, to the equation 


o o 

-^po{k,T) = -Xk-^poik,T) - |fc|“J'[|a:r%(a;,r)]. 


(18) 


Eq.dnj cannot be exactly solved for arbitrary x] we look for a solution in the limit of small |fc|. It can be represented 
by a stable distribution but only in a limited range of 0 [s^. The solution valid for an arbitrary large 0 contains a 
power-law tail but its shape at small |a:| is different from that for the stable distribution [33 |. In this paper, we assume 
that the solution has a scaling form that, in the asymptotic regime, is defined by the parameter a, 


Po{x,t) = a{T)po{a{T)x) = a{T) “|x| “ \ 


(19) 


and the validity of a tentative assumption (1191) will be verified and the unknown function a(r) determined by inserting 
Eg. lfT^ to Eg. dTBl) . For that purpose, both Fourier transforms in Eo. dTSl) have to be evaluated in a limit of small |A:|. 
Transform in the last term is given by the integral (nni), 


a(r) 


|a;| ^Po{a{T)x)dx = a(T)® 


|a;| ^po(x)dx = airfho, 


( 20 ) 


where hg =const. By the evaluation of the Fourier transform from po{x^ r) we apply the same procedure as in Sec.II 
and obtain 


Po{k,T) = 1 - 


a(T) “tt 


r(a -I- 1) sin(a;7r/2) 

Inserting the above expression to Eo. lfT^ produces a differential equation. 


\k\^ 


dM — Aa H- air) 

ac 


a+0+l 


= 0 , 


which can be solved by separation of the variables: 


a(T) = A 




-l/(a+9) 


( 21 ) 


( 22 ) 


(23) 


in this equation, A = and c is defined by Eq. (US. The final asymptotic solution is given by Eq. dUD 

and this expression is exact in respect to time- and position-dependence; the constant Hq contains details of the 
solution for all x and cannot be determined by the above procedure. The case 0 = 0 is exceptional: then Eq. (HU can 
be solved exactly for arbitrary x and the solution resolves itself to the stable distribution Q. 


III. RELAXATION TO STATIONARY STATES 

Stochastic properties of the systems for which particle is subjected to the Gaussian white noise and remains under 
the influence of an external deterministic force are usually described by the Fokker-Planck equation. The expansion 
of its solution to the eigenfunctions of the Fokker-Planck operator contains single modes which decay exponentially in 
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time to a stationary state representing a Gibbs-Boltzmann distribution and the temperature is related to the diffusion 
coefficient by the fluctuation-dissipation theorem [T^ . Generalisation of that equation takes into account the memory 
effects by introducing a fractional derivative over time (the fractional Fokker-Planck equation). In this generalised 
description, a different relaxation pattern of the single modes emerges: the decay is governed by the Mittag-Leffler 
function [l2j| and, in the case of the harmonically bound particle, also the variance d ecay s according to that pattern 
[ 3 ^. The Mittag-Leffler relaxation naturally emerges in the subordination technique [T^ . 

The system, relaxation properties of which are analysed in the present paper, is defined by the stochastic equation 
(|51) . It includes the random force in a form of the Levy stable distribution and the non-Markovian nature of the 
process is taken into account by applying a subordination technique: the process - described in Sec.II and given by 
Po{x,t) - is subordinated to the random time. This time is defined by a one-sided, maximally asymmetric stable 
Levy distribution h'{T,t) = Lpir), where 0 < /3 < 1, which vanishes for r < 0. Since the processes x{t) and r(t) are 
statistically independent, the distribution as a function of the physical time is given by the total probability formula. 


p{x,t)= po{x,T)h{T,t)dT, (24) 

^0 

where h{T^t) = denotes an inverse subordinator El- In the next subsection, we derive that distribution 

for the linear case. Note that stationary state is completely determined by the distribution pq(x,t), i.e. by the noise 
type and the potential, whereas the subordinator, as a source of the subdiffusion, influences only the relaxation to 
that state. 


A. Linear case 


Let us consider the case 7 = 1 -I- d which condition means that the effective deterministic force in the first equation 
dU is linear and the subordinated process is governed by Ea. lflTll . Note that then the force Fd in Eq.® may be 
nonlinear. The density p{x,t) follows from Eg. (1241) and can be evaluated by a direct integration. For that purpose, 
we insert the expression (1211) to Ea. (l24l) and calculate the Fourier and Laplace transform taking into account that the 
Laplace transform from ft,(r, t) is given by h{T, u) = exp(— tu^). This procedure yields 


p{k,u) 


1 

-/ a(T)' 

u Jo 


and the expression for a(T) “ results from Eg. (1231) by applying a direct exponentiation: 


(25) 


°° /_1 \j 

A“a(r)-“ = 1 + 51 




X{a+e)jT 


J'- 




j=l •' i=0 

where q = a/{a + 9). The integration yields the transformed density, 

0-1 


p{k, u) = — — cA 
u 


1 , i-iy 

u ^ + A(a + e)j AJp 




(26) 


(27) 


The above transform can easily be inverted when we realise that the u—dependent term corresponds to a Mittag-Leffler 
function that is defined by the series 


E0{x) 


00 


E 

n—O 


r(l+/3n) 


(28) 


and resolves itself to an ordinary exponential for /3 = 1. The final expression for the asymptotic form of the density 
in the Eourier space reads. 


where a function 


p{k,t) = l-cA “[1 - x(t)]|fc|“, 


xit) 


00 / ly- J -1 

E —^E0[-X{a + 9)jt>^] n(9 - *)> 

j=l i=0 


(29) 


(30) 
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FIG. 1: (Colour online) Relaxation function x{i) calculated from Eg. (1301) for a 


1.5, p = 0.5, A = 1 and a few values of 6. 


involving a series of the Mittag-Leffler functions, determines the relaxation to the stationary state and x(oo) = 0. 
Inversion of the Fourier transform yields a time-dependence of the asymptotic solution, 

p(x,t) = Al-“[l-x(f)]N-“-i (|x|»l). (31) 

Notice that x(0 does not depend on ho and is exact, in contrast to the intensity of p(x^ t)\ the latter quantity depends 
on the shape of the density at small \x\ and is not determined by the present method. The long-time behaviour of 
the relaxation function x(i) agrees with that of the Mittag-Leffler function and then x(0 but the effective 

relaxation time may strongly depend on the parameter 9. Fig.l presents x{t) calculated from Eg. ( 1501 ) for a few values 
of 9. The Mittag-Leffler function was determined from Eg. (1^51) for small t and from the expansion 


Ei3{x) 


h r(i - fin) 


(32) 


for large t; the convergence could not be achieved for intermediate values [s^ . Despite the same asymptotic form, the 
relaxation time substantially differs for various 9: it is very large for 9 < 0. 


The time-dependent density distribution for a stochastic dynamical system satisfies a Fokker-Planck equation which, 
in the presence of jumps and/or long waiting times, is fractional. This equation can be derived from the generalised 
master equation and CTRW [l^; it has been applied, in particular, to describe chaotic systems characterised by a 
hierarchical structure of the trajectories [s^. The equation in the latter case possesses, in general, a variable diffusion 
coefficient. The fractional Fokker-Planck equation has a stochastic representation in a form of equations similar to 
O- This has been proved for the homogeneous case with an arbitrary potential and that equivalence can be 
generalised to the nonhomogeneous case. Indeed, we demonstrate in the Appendix A that the function p(x, t), Eg. (1311) . 
satisfies the following Fokker-Planck equation, 


dpjx, t) 
dt 


= oA 


d 


do 


A^[a;p(a:,t')] + p{x,t)) 


d\a 


where the Riemann-Liouville fractional derivative is defined by the integral operator, 


(33) 




1 d 


dt' 


fit') 


TiP)dtJo 


(34) 


B. Nonlinear case: a numerical analysis 

In this subsection, we consider the general case when the effective deterministic force in the first equation (|3]) is 
a nonlinear function of a: (7 1 -I- 9). To evaluate the density p{x,t) from Ea. (l24ll we need the distribution for the 
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FIG. 2: (Colour on line) The time-evolution of the density distribution for the subordinated process corresponding to two 
values of 7 and calculated for a = 1.5 and 9—1 from Eq.®. The curves in the upper part correspond to r =0.05, 0.3, 1, 5 
and 20 (from bottom to top). The straight red lines mark the dependence and x~°‘ Averaging was performed over 

10 ^ trajectories. 


subordinated process, pq{x,t), i.e. the solution of Eq.([7]). This distribution cannot be found analytically, even for 
0 = 0, if 7 = 7 — 07 ^ 1 . Therefore, we must resort to a numerical analysis details of which are presented in Appendix 


B. 


First, we consider Eq. to get a qualitative picture of the solutions for various values of 7 . Just a simple analysis 
of Eq.® indicates that the cases 7 > 1 and 7 < 1 must be different. The particle, initially positioned at a; = 0, is 
subjected only to the random force at the beginning of the time-evolution since the first term in Eq.® can be neglected 
for 7 > 0. Hence the density distribution corresponds to the case without any potential and pq{x,t) oc 
However, the density approaches the stationary limit, which is of the form differently for both cases: the 

tail becomes steeper (7 > 1) or flatter (7 < 1) during the evolution. The density distributions for those cases are 
presented, as a function of time, in Fig.2; we assume, from now on, A = 1. If 7 is large, the asymptotic distribution, 
|a;|-a- 7 ^ is visible already at a very short time but restricted to a large distance. The stationary distribution is 
reached at r = 1 and the curves for the intermediate cases show a widening of the stationary tail. In contrast to this 
picture, the case 7 = 0.3 exhibits a gradual change of the slope for the entire tail and the stationary distribution is 
reached at r = 10 . 


We consider both cases separately and assume, for the large 7 , the finite variance: a -|- 7 > 3. Then 

p{k, t) Ri 1 — a{t)k^ = 1 — aek^{l — x(t)), (35) 

where a{t) is a function describing a time-dependence of the asymptotic density which converges to its stationary 
value tte = a(oo); on the other hand, we have a(0) = 0 for the initial condition in the form of the delta function. The 
relaxation function 


X{t) = 1 - a{t)/ae = 1 - {x'^){t)/{x'^){oo) 


( 36 ) 










FIG. 3: (Colour on line) The relaxation function for a = 1.9, /3 = 0.5, 7 = 3 and three values of 6. The solid lines on the right 
side mark an exponential, with = 0.33, 0.4 and 0.14 (from left to right). The solid line on the left side, fitted to 

the data for 9 = 1, marks a function Stars correspond to the case 7 = 4 and 9=1. Averaging was performed over 10"^ 

trajectories. 


will serve as an estimator of the convergence rate to the stationary state. It is shown in Fig.3 for a few values of 
9. In contrast to the linear case, the relaxation is not smooth: it exhibits two regions of a very different slope. At 
a short time, x(f) diminishes relatively slowly whereas at large t it rapidly falls according to the exponential rate. 
This picture emerges for both positive and negative 6 but the break down of the curve is especially pronounced for 
0 = —0.5. The case 9 = 0 also indicates the exponential fall but the difference between the behaviour for small and 
large t is less striking. The decline at small t and for d ^ 0 obeys, approximately, a power-law dependence; the slope 
is more gentle than for the linear case. We conclude that instead of the power-law pattern, present for the linear case, 
the exponential relaxation is observed at long times. The curves in Fig.3 correspond to different effective potentials 
in Eq.([3|) and one can ask whether this is a reason of their different shapes. To check that, we performed a calculation 
for 7 = 4 and 6 = 1 which corresponds to the same effective potential as the case 7 = 3 and 0 = 0. Results, presented 
in Fig.3, indicate a clear disagreement of those cases and the relaxation pattern is determined by 9: the curves for 
the parameter sets 7 = 4, 0 = 1 and 7 = 3, 0 = 1 actually coincide. 


Estimation of the decay rate to the stationary state in terms of the function x(t), Ea. (l5Sl) . is not possible for 7 < 1 
(7 < 0 + 1) since then the variance is infinite. Instead, we consider the time-evolution of the slope of the asymptotic 
distribution, Evaluation of the tails is numerically difficult - in respect to the precision and the calculation time 
- if one directly solves Eq.(|31). On the other hand, the distributions in the operational time for the process defined by 
Eq.(|3]) are easier to handle and the tails can be precisely determined. Therefore, we calculated from Eq.(|3]). As 
the first step, the slope and the relative intensity of Cmir), were evaluated from Eq.(|ni); the densities 

are presented in the upper part of Eig.2. Knowing those quantities, the final asymptotic distribution p{x,t) can be 
obtained by a direct integration by means of Eg. (Oil) . 

POO 

p{x,t)= / (37) 

Jo 


and we restrict our analysis to the case (3=1/2 for which the inverse subordinator has a simple form. 
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FIG. 4: Slope of the density po{x,t) calculated from Eq.® for a = 1.5, /3 = 0.5, 7 = 0.3 and 6 — 1 (upper part). The 
lower part presents the slope as a function of the physical time, calculated from Eq. (f37l). Lines mark the functions r ° and 
^-0.049 upper and lower part, respectively. Each point was obtained from the density evaluated by averaging over 10^ 

trajectories. 


Slopes of the asymptotic distributions characterising both processes are presented in Fig.4. They have different shapes 
as functions of time: Prn{T) is almost constant at small times and declines as a power-law reaching its stationary value 
at T = 20. The integration (IWl) brings about a smoothing effect: /i(<) has a power-law shape for all times and the 
index /i is smaller than for the subordinated process at a corresponding time; the stationary state is reached at t = 10. 


IV. SUMMARY AND CONCLUSIONS 

We have considered a dynamical process defined by the stochastic stimulation - governed by the Levy stable 
statistics - and the external potential in the form of both harmonic and strongly non-linear oscillator. It has been 
assumed that the particle may be trapped for a long time and the trapping time varies with the position, as a result of 
a nonhomogeneous medium structure. The problem is formulated in terms of the subordination technique where the 
intensity of the random time in the subordinator is position-dependent. The asymptotics of the stationary state has 
been derived for the case of the non-linear deterministic force. It has a power-law form; the slope depends on both the 
stability index of the driving noise and the potential but is independent of the memory characteristics: /3 and 9. The 
relaxation process to that stationary state has been analysed for both harmonic and non-linear case. For the linear 
case (the effective potential is quadratic), the decay rate is observed at long times but the value of the relaxation 
function x strongly varies as a function of 0; long trapping times near the origin (6* < 0) correspond to a long time 
of approaching the stationary state. In contrast to the case 9 = 0, where relaxation obey the Mittag-Leffler pattern, 
x(t) for the general case involves a linear combination of the Mittag-Leffler functions. The density distribution for 
the linear case satisfies the integro-differential equation with a complicated kernel that resolves itself to a fractional 
equation with a variable diffusion coefficient in a limit of long time. 

The density distributions for a nonlinear effective force exhibit a complicated and differentiated structure which 
qualitatively depends on the potential slope. The strong potentials mean a rapid development of the stationary 
slope whereas the weak potentials are characterised by a gradual flattening of the entire tail. In the former case, 
the relaxation function is always exponential for long times: the Mittag-Leffler pattern breaks down and the Debye 
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relaxation is observed. At smaller times, the decay is slow and weakly depends on 9. Numerical calculations indicate 
that the relaxation pattern is sensitive rather on the spacial nonhomogeneity of the memory, parametrised by 9, than 
on a slope of the effective deterministic force. 

The results indicate that the relaxation for a complex system ~ in which long tails of the random force distribution 
coexist with a position-dependent trap distribution - to a stationary state depends on both an external potential 
and a structure of the trap positions. When the non-exponential pattern is observed, as it is the case for the linear 
systems, the latter characteristics may strongly influence the relaxation time. However, the linearity property refers 
not to the external potential alone but results from an interplay between this potential and the trap distribution. 
Therefore, though the nonlinearity changes the relaxation pattern to the exponential, presence of the Mittag-Lefller 
pattern in the experimental data does not necessarily imply that the linear deterministic force would be correct in a 
dynamical description of the case of the uniform trap distribution. Finally, we emphasise the non-homogeneous trap 
distribution can make the variance of the process value finite. The divergent variance for the Levy flights is often 
unphysical and we demonstrated that this difficulty may not emerge even for relatively weak potentials if 9 is negative 
and sufficiently small. 


APPENDIX A 


In the Appendix A, we demonstrate that the density p{x,t), Eg. (1311) . satisfies the Fokker-Planck equation (l33t . 
Let us consider the following integro-differential equation. 




X^[xp{x,t')] + -^{\x\ %{xX)) 


dt' 


(Al) 


and take the Fourier-Laplace transform which we rewrite as Fi = K{u){F 2 + F^). Next, we apply Ea. (l27ll and obtain 


Fi ^ C[—p{k,t)] =-cA 




Similarly, the first term on rhs reads 


F 2 =—Xk—p{k,u) = cXaA 
ok 


“ j\ u°-\-X{a-\-0)j 


i + u/3-i ni=d(9-*) 

n, 


(A2) 


j! u9 + X{a + 9)j 


(A3) 


and the last one, 

nOO 

F 3 = -\k\‘^F[\x\-‘^p{x,u)] = = -\k\^hoA^ 

Jo 

where q' = —9/(a + 9) = q — 1. Denoting the sum in Ea. (IA2l) by Si, 

i-^y niZoid-y 


^ Al + ^)j 






j\ u9 + A(a -I- 9)j 


and, similarly, the sum in Ea. (IA3l) by S 2 , we get 

F2 +F3 = cXaA-^\k\^u^-\Si - 82)-, 

in the above equation we took into account that A® = ^^A~°‘. The expression Si — S 2 can be simplified by changing 
the product index in the second sum: 


(A4) 

(A5) 

(A 6 ) 




(-1V 1 (-IP niZoi^-i) 


—^ ?! X(a0)j 




^ U - 1)! + X{a -f 9)j 


: = q-^S. (A7) 


Dividing Fi by F 2 + F 3 yields the expression for the Laplace transform from the kernel, 

K(u) = _i_i±^^l -/5 

^ ’ X{a + 9) S ’ 


(A 8 ) 
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which can be evaluated by using the identity, 


-*) = -!• (A9) 

1=1 i =0 

Finally, K{u) = , inversion of this transform produces the operator (IMl) and Ea. dAll) takes a form of the fractional 

equation (1551) . 


APPENDIX B 

In the Appendix B, we present numerical details of the numerical evaluation of the stochastic trajectories from 

Eq.®. 

The time-evolution of a stochastic system which is represented by the subordination equations can be numerically 
evaluated, namely the process value as a function of the physical time can be found, without an explicit inversion of 
the subordinator [^, . We solve a coupled set of equations ® the first of which is an ordinary Langevin equation 

and can be discretized, 


x{t„) = a;(r„_i) -t Fd[a;(T„_i)]Ar -b AyyAr^/". (Bl) 

The second equation ®, in turn, represents a strictly increasing stable motion (since g{x) > 0), defined as 

r(t) = inf{r : t^r) > t}, (B2) 


and the discretized equation reads 

tiTn) = -|-5(x(t„_i))A^At^/^. (B3) 

Increments of the stable distributions, Ar] and A^, are sampled by means of the standard procedure [dlj. Operationally, 
x{T) follows from the integration of both equations as long as t > T; finally, the definition (IB2I) is applied. In the 
calculations presented in Eig.2-4, the step Ar = 10“'* was used. 
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